#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*====================================================*/

#define DEBUG 1
#define MAXX 64

float hprint(unsigned int  *hist, int nbin, char *title) {
  int i,j,k,v,s,s0,s1,s2,s3;    
  int ybin,y5,miny=15,maxy=25,mod,entry;
  float sum,sum2,mean,mean2,ey2,rms,x;
  char st[MAXX+1];
  /*----------------------------------------------------*/
  int sum_bin,max,max1;
  unsigned int pd[MAXX];
  printf("\n");

  /*---  X scale (rebin) ---*/

  for (i=1;i<32768;i++)  if ((nbin/i)<=MAXX) break;
  sum_bin=i;
  if (DEBUG) printf(" Xscale: nbin=%d sum_bin=%d \n",nbin,sum_bin);

  /*---  Y scale ---*/

  i=0; j=0; max=0; max1=0; entry=0;
  while(1) {
    pd[i]=0;
    for (k=0;k<sum_bin;k++) { 
      if (j<nbin) { pd[i]+=hist[j]; if (hist[j]>max1) max1=hist[j]; } else pd[i]+=0; 
      j++; 
    }
    if (pd[i]>max) max=pd[i];  entry+=pd[i];
    if (DEBUG>1) 
      printf(" pd[%d]=%d max=%d \n",i,pd[i],max);
    i++; if (i==MAXX) break;
  }

  for (j=10;j>0;j--) {
    mod=j%10;
    for (i=maxy;i>miny;i--) {
      ybin=max/i+1; 
      if (DEBUG>1) 
        printf(" Yscale: maxY=%d, max=%d, ybin=%d mod=%d\n",i,max,ybin,mod);
      if ((ybin%10)==mod) goto a;
    }
  }
a:
  if      (max<=25) { maxy=max; ybin=1; }
  else if (i!=miny) { maxy=i; 
                      if (mod) {   if (ybin>=7)  ybin+=10-mod; 
			maxy=max/ybin+1;
			if (maxy<16) maxy=16;
		      } 
                    } 
  else              { ybin=max/maxy+1; }
  if (ybin==0)  ybin=1;

  /*-------------  Statistics   ------------------------*/
  mean=0; rms=0;
  if (entry>0) {
    sum=0; sum2=0;
    for (i=0;i<nbin;i++) {     x=i+0.5;
      sum +=(float)hist[i]*x;
      sum2+=(float)hist[i]*x*x;
    } 
    mean=sum/entry; mean2=mean*mean; ey2=sum2/entry;
    if (DEBUG>1) 
      printf(" sum=%f, sum2=%f, \n mean=%f, mean2=%f, E(y**2)=%f\n"
	     ,sum,sum2,mean,mean2,ey2);
    rms=(float)sqrt((double)((float)nbin/(float)(nbin-1)*(ey2-mean2)));
  }
  printf(" Title: %s\n",title);
  printf(" Entry=%d, Max=%d(%d), Mean=%.2f, RMS=%.2f \n"
          ,entry,max,max1,mean,rms);
  /*----------------------------------------------------*/
  if (DEBUG>1) printf(" Yscale: maxY=%d  max=%d  ybin=%d\n",maxy,max,ybin);
  max=maxy*ybin; y5=ybin/5;
  if (DEBUG) 
    printf(" Yscale: maxY=%d  max=%d  ybin=%d  y5=%d \n",maxy,max,ybin,y5);

  /*----- boba -----------------------------------------*/

  for (j=1;j<=maxy;j++) {
    s=max-ybin*j;   st[MAXX]=0;
    s0=s+y5; s1=s0+y5;  s2=s1+y5;  s3=s2+y5;
    for (i=0;i<MAXX;i++) {   /*  fill line  */
      v=*(pd+i);      st[i]=' ';
      if (0<v && v<=ybin && j==maxy && ybin!=1)   {st[i]='.';  }
      if (v>s)       {st[i]='.';     }
      if (v>s0)      {st[i]=',';     }
      if (v>s1)      {st[i]='+';     }
      if (v>s2)      {st[i]=';';     }
      if (v>s3)      {st[i]='X';     }
      if (DEBUG>2) 
	printf(" v=%d, s0123=%d %d %d %d %d ybin=%d, y5=%d, i=%d st[i]=%c \n"
	      ,v,s,s0,s1,s2,s3,ybin,y5,i,st[i]);
    }
    printf("%8d   I%s\n",s+ybin,st);
  }
  for (i=0;i<MAXX;i++) { st[i]='0'+i%10; }
  printf("            %s\n",st);
  printf("   ");
  for (i=0;i<MAXX/9;i++) { printf("%10d",10*i*sum_bin); }
  printf("\n");
  return mean;
}
/*=============================================================*/
/*                                                             */
/*=============================================================*/
float gauss()
{
  int i;
  float s = -6.0;
  for(i=0;i<12;i++) s += ((float)rand())/((float)RAND_MAX);
  return s;
}

void mfill(unsigned int *hist, int n, float x, float w) {
  int ich;
  ich=(int)x; 
  if (ich>=n) { ich=n-1;/* printf(" mfill: Over:  n=%d x=%f \n",n,x);*/}
  if (ich<0)  { ich=0; /*  printf(" mfill: Under: n=%d x=%f \n",n,x);*/}
  hist[ich]+=(int)w;
  if (DEBUG>2) printf(" MFILL: x=%f, ich=%d w=%d \n",x,ich,w);
  /* printf(" MFILL: hist=%p Len=%d \n",hist,n); */
}

void mreset(unsigned int *hist, int n) {
  int i;
//  printf("mreset:: hist=%p len=%d \n",hist,n);
  for (i=0;i<n;i++)  hist[i]=0;
}
void mcopy(unsigned int *hist1, unsigned int *hist2, int n) {
  int i;
  for (i=0;i<n;i++)  hist2[i]=hist1[i];
}

/*=============================================================*/
/*                                                             */
/*=============================================================*/
float hstati(int sw, unsigned int  *hist, int nbin) {
  int i,entry;
  float sum,sum2,mean,mean2,ey2,rms,x;
  /*-------------  Statistics   ------------------------*/
  entry=0; sum=0; sum2=0;
  for (i=0;i<nbin;i++) {     x=i+0.5;
    entry+=hist[i];
    sum +=(float)hist[i]*x;
    sum2+=(float)hist[i]*x*x;
  } 
  if (entry==0) return 0;
  mean=sum/entry; mean2=mean*mean; ey2=sum2/entry;
  if (DEBUG>1) 
    printf(" sum=%f, sum2=%f, \n mean=%f, mean2=%f, E(y**2)=%f\n"
	   ,sum,sum2,mean,mean2,ey2);
  rms=(float)sqrt((double)((float)nbin/(float)(nbin-1)*(ey2-mean2)));
  /*printf(" Entry=%d, Mean=%.2f, RMS=%.2f \n",entry,mean,rms);*/  
  /*----------------------------------------------------*/
  if     (sw==1)  return (float)entry; 
  else if (sw==3) return rms;   
  else            return mean;  

}


